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Abstract 

We consider the problem of one dimensional (ID) phase retrieval, namely, recovery of a ID signal from the 
magnitude of its Fourier transform. This problem is ill-posed since the Fourier phase information is lost. Therefore, 
prior information on the signal is needed in order to recover it. In this work we consider the case in which the 
prior information on the signal is that it is sparse, i.e., it consists of a small number of nonzero elements. We 
propose a fast local search method for recovering a sparse ID signal from measurements of its Fourier transform 
magnitude. Our algorithm does not require matrix lifting, unlike previous approaches, and therefore is potentially 
suitable for large scale problems such as images. Simulation results indicate that the proposed algorithm is fast and 
more accurate than existing techniques. 

I. Introduction 

Recovery of a signal from the magnitude of its Fourier transform, also known as phase retrieval, is of great interest 
in applications such as optical imaging lUl, crystallography [2J, and more [3|. Due to the loss of Fourier phase 
information, the problem (in ID) is generally ill-posed. A common approach to overcome this ill-posedeness is to 
exploit prior information on the signal. A variety of methods have been developed that use such prior information, 
which may be the signal's support (region in which the signal is nonzero), non-negativity, or the signal's magnitude 
(H, im. A popular class of algorithms is based on the use of alternate projections between the different constraints. 
In order to increase the probability of correct recovery, these methods require the prior information to be very 
precise, for example, exact/or "almost" exact knowledge of the support set. Since the projections are generally 
not onto convex sets, convergence to a correct recovery is not guaranteed [6|. A more recent approach is to use 
matrix-lifting of the problem which allows to recast phase retrieval as a semi-definite programming (SDP) problem 
lITl. The algorithm developed in fT\ does not require prior information about the signal but instead uses multiple 
signal measurements (e.g., using different illumination settings, in an optical setup). 

In order to obtain more robust recovery without requiring multiple measurements, we develop a method that 
exploits signal sparsity. Existing approaches aimed at recovering sparse signals from their Fourier magnitude belong 
to two main categories: SDP-based techniques llSl.llQl. lfTOi . lfTTl and algorithms that use alternate projections (Fienup- 
type methods) [12|. Phase retrieval of sparse signals can be viewed as a special case of the more general quadratic 
compressed sensing (QCS) problem considered in |8|. Specitically, QCS treats recovery of sparse vectors from 
quadratic measurements of the form yi = x'^AiX, i = 1, . . . ,N, where x is the unknown sparse vector to be 
recovered, are the measurements, and Aj are known matrices. In (discrete) phase retrieval, A = FJFi where 
Fj is the zth row of the discrete Fourier transform (DFT) matrix. QCS is encountered, for example, when imaging 
a sparse object using partially spatially-incoherent illumination [8]. 

A general approach to QCS was developed in lH based on matrix lifting. More specifically, the quadratic 
constraints where lifted to a higher dimension by defining a matrix variable X = xx^. The problem was then 
recast as an SDP involving minimization of the rank of the lifted matrix subject to the recovery constraints as 
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well as row sparsity constraints on X. An iterative thresholding algorithm based on a sequence of SDPs was then 
proposed to recover a sparse solution. Similar SDP-based ideas were recently used in the context of phase retrieval 
I0,|[IO|. However, due to the increase in dimension created by the matrix lifting procedure, the SDP approach is 
not suitable for large-scale problems. 

Another approach for phase retrieval of sparse signals is adding a sparsity constraint to the well-known iterative 
error reduction algorithm of Fienup (TT\ . In general, Fienup-type approaches are known to suffer from convergence 
issues and often do not lead to correct recovery especially in ID problems; simulation results show that even with 
the additional information that the input is sparse, convergence is still problematic and the algorithm often recovers 
erroneous solutions. 

In this paper we propose an efficient method for phase retrieval which also leads to good recovery performance. 
Our algorithm is based on a fast 2-opt local search method (see 1 13] for an excellent introduction to such techniques) 
applied to a sparsity constrained non-linear optimization formulation of the problem. We refer to our algorithm 
as GESPAR: GrEedy Sparse PhAse Retrieval. Sparsity constrained nonlinear optimization problems have been 
considered recently in [14J; the method derived in this paper is motivated - although different in many aspects - 
by the local search-type techniques of L14|. In essence, GESPAR is a local-search method, where the support of the 



sought signal is updated iteratively, according to selection rules described in detail in Section III A local minimum 
of the objective function is then found given the current support using the damped Gauss Newton algorithm. We 
demonstrate through numerical simulations that the proposed algorithm is both efficient and more accurate than 
current techniques. 



The remainder of the paper is organized as follows. We formulate the problem in Section |lT] Section III describes 



our proposed algorithm in detail. Numerical performance is illustrated in Section IV 



II. Problem Formulation 
We are given a vector of measurements y E M^, that corresponds to the magnitude of an point DFT of a 



vector X G M^, i.e.: 



yi 



2,rj(m-l)(i-l) 

Xme " 



l = l,...,N, (1) 



m=l 

where x is constructed by zero padding of a vector x G with elements Xj, i = 1, 2, . . . , n. In the simulations 
section we consider the setting N = 2n which corresponds to oversampling the DFT of x by a factor of 2. 
More generally, we assume that N > 2n — 1. This allows to determine the correlation sequence of x from the 
given measurements, as we elaborate on further below. Denoting by F G £,NxN jj^g DFT matrix with elements 
exp I — ^iMn^LZpH—ll'j ^ we can express y as y = |Fx|, where | • | denotes the element- wise absolute value. The 
vector X is known to be s-sparse, that is, it contains at most s nonzero elements, and in addition its last n elements 
are zero. Our goal is to recover x given the measurements y and the sparsity level s. 

The mathematical formulation of the problem that we consider consists of minimizing the sum of squared errors 
subject to the sparsity constraint: 

minx E^Ii(|F.x|2-y2)2 

s.t. ||x||o < s, ^2) 
supp(x) C {1,2, ...,n}, 
X G M^, 

where Fj is the ith row of the DFT matrix F, and || • ||o stands for the zero-"norm", that is, the number of nonzero 
elements. Note that the unknown vector x can only be found up to trivial degeneracies that are the result of the 
loss of Fourier phase information: circular shift, global phase, and signal "mirroring". 

To aid in solving the phase retrieval problem we will rely on the fact that the correlation sequence of x (the first 
n components of x) can be determined from y. Specifically, let gm = Y17=i XiXi+m,rn = — (n — 1), . . . ,n — 1 
denote the correlation sequence. Note that {gm} is a sequence of length 2n — 1. Since the DFT length satisfies 
> 2n — 1, we can obtain {gm} by the inverse DFT of the squared Fourier magnitude y. Throughout the paper, 
we assume that no support cancelations occur in {gm}, namely, if Xj 7^ and Xj ^ for some i, j, then g\i^j\ 7^ 0. 
When the values of x are random, this is true with probability 1. This fact is used in GESPAR in order to obtain 
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initial information on the support of x. 

The information on the support is used to derive two sets, Ji and J2 from the correlation sequence {g-m} in the 
following manner. We denote by Ji the set of indices known in advance to be in the support, from the autocorrelation 
sequence. To derive the set Ji, note that due to the existing degree of freedom relating to shift-invariance of x, 
the index 1 can be assumed to be in the support, thereby removing this degree of freedom; as a consequence, the 
index corresponding to the last nonzero element in the autocorrelation sequence is also in the support, i.e. 

imax = 1 + argmax{i : / 0}. 

i 

Therefore, Ji = {I, imax}- 

Next, we denote by J2 the set of indices that are candidates for being in the support, meaning the indices that 
are not known in advance to be in the off-support (the complement of the support). Specifically, J2 contains the set 
of all indices k G {1, 2, . . . , n} such that gj^^i / 0. Obviously, since we assume that = for A; > n, we have 
-^2 ^ {1, 2, . . . , n}. As a concrete example, consider the signal x = (2, 0, 0, —1, 0, 1.5)"'^. The corresponding 11 point 
autocorrelation function gm is then given by gm = (—3, 0, —2, 1.5, 0, 7.25, 0, 1.5, —2, 0, —3)-^. The set Ji is therefore 
Ji = {1, 6}. Next, by examining the zeros of g^, and using our assumption of no support-cancelations, we deduce 
that there are no two non-zero elements in x 7^ and Xj 7^ such that \i — j\ = 1, 4. Therefore, forcing the first 
element in x to be non-zero, which removes the shift-invariance degeneracy, immediately implies that X2 = = 0. 
In this way J2 is determined as J2 = {1,3,4,6}. Defining Ai = 3ft(Fi)^3f?(Fi) + 9(Fi)'^9(Fi) G M^^^ and 
Ci = yf for i = 1, 2, . . . , A^, problem ^ along with the support information can be written as 

minx /(x) = Z]iIi(x^AiX - a)^ 
s.t. ||x||o < s, 

Ji C supp(x) C J2, 

X G M^, 



(3) 



which will be the formulation to be studied. 

In the next section, we propose GESPAR - an iterative local-search based algorithm for solving Q. We note 
that although in the context of phase retrieval the parameters Aj, Ji, J2 have special properties (e.g., A^ is positive 
semidefinite of at most rank 2, | Ji| = 2), we will not use these properties in GESPAR. Therefore, our approach is 
capable of handling general instances of Q with the sole assumption that Aj is symmetric for any i = 1,2, . . . , N. 

III. Greedy Sparse PhAse Retrieval (GESPAR) 

A. The Damped Gauss-Newton Method 

Before describing our algorithm, we begin by presenting the damped Gauss-Newton (DGN) method fT5l,fT6l that 
is in fact the core step of our approach. The DGN method is invoked in order to solve the problem of minimizing 
the objective function / over a given support 5 C {1, 2, . . . , n} {\S\ = s): 

min{/(U5z) : z G R^}, (4) 

where U5 G M"^* is the matrix consisting of the columns of the identity matrix I^v corresponding to the index set 
S. With this notation, Q can be explicitly written as 

min |5(z) = J^(z^U^AiU5Z - af : z G | . (5) 

Problem (|5]l is a nonlinear least-squares problem. A natural approach for tackling it is via the DGN iterations. 
This algorithm begins with an arbitrary vector zq. In our simulations, we choose it to be an uncorrelated random 
Gaussian vector with zero mean and unit variance. At each iteration, all the terms inside the squares in g(z,) are 
linearized around the previous guess. Namely, we write g{'L) from ([S]) as: 

N 

g{z) = ^h^{z) (6) 

i=l 
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Algorithm 1 DGN for solving ^ 
Input: {Ai,Ci,S,e,L). 

Aj G M^^^, i = 1, 2, . . . , - symmetric matrices. 

CieR,i = l,2,...,N. 

5 C {1, 2, . . . , n} - index set. 

e - stopping criteria parameter. 

L - maximum allowed iterations. 

Output: z - an optimal (or suboptimal) solution of Q. 
Initialization: Set = \J'gAiUs,to = 0.5, zq a random vector. 

General Step k{k > 1): Given the iterate z^-i, the next iterate is determined as follows: 

1. Gauss-Newton Direction: Let y^. be the solution of the linear least squares problem [8j given by: 

Yk = (M^M)-iM^b 

with the ith row of M being Mj = 2(BiZfc_i)-^, and bi = Ci + z|l^BjZfc_i for i = 1, 2, . . . , A^. The Gauss-Newton 
direction is = - Zk-i. 

2. Stepsize Selection via Backtracking: set u = min{2tfc_i, 1}. Choose a stepsize tk as tj. = {^)"^u, where m is 
the minimal nonnegative integer for which g{zk-i + (^)™tidfe) < g{zk^i). 

3. Update: set z^ = Zjt_i + t^d^. 

4. Stopping rule: STOP if either \\zk — Zk-i\\ < e or k > L. 



, with hi{z) = z^BjZ — Ci, where Bj = U^AjUs. At each step we replace hi by its linear approximation around 

hi /ij(zfc_i) + V/ii(zfc_i)^(y - Zfc_i) 
= z^_iBjZfc_i - Q + 2(BjZfc_i)^(y - Zk-i) 
We then choose y^ to be the solution of the problem 

N 

min V(zi'_iB,Zfe_i - a + 2(B,Zfe_i)'^(y - Zk-i)f. (8) 
y 

i—i 

Problem ([8]l can be written as a linear least squares problem 

yfc = argmin ||My - b||^ (9) 

with the ith row of M being Mj = 2(BjZfc_i)^, and bi = Ci + z|l^BjZfc_i for i = 1, 2, . . . , A^. The solution y^ 
is given by M, i.e. yk = (M^M)~^M^b. We then define a direction vector as d^ = y^ — Zfc_i. This direction is 
used to update the solution with an appropriate stepsize designed to guarantee the convergence of the method to a 
stationary point of g{z). The stepsize is chosen via a simple backtracking procedure. Note that the method requires 
M to be of a full column rank. In our numerical experiments we encountered cases in which this assumption was 
not valid, and in these cases, our implementation chose one of the optimal solutions of the corresponding least 
squares problem. We noticed that these cases had negligible effect on the results 

Algorithm [T] describes the DGN method in detail. In our implementation the stopping parameters were chosen 
as e = 10"^ and L = 100. 



B. The 2-opt Local Search Method 

The GESPAR method consists of repeatedly invoking a local-search method on an initial random support set. 
In this section we describe the local search procedure. At the beginning, the support is chosen to be a set of s 
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Algorithm 2 2-opt 
Input: (Ai,Q, Af,P). 

Aj G M^^^, i = 1, 2, . . . , - symmetric matrices. 
Ci G]R,i = 1,2,..., iV. 
M, P - positive integers. 

Output: X - a suggested solution for problem Q. 
T - total number of required swaps. 



1) Initialization: 

a) Set T = 0. 

b) Generate a random index set S'o(|S'o| = s) satisfying the support constraints (Ji C 5o ^ J2). 

c) Invoke the DGN method with parameters (Aj, Cj, /So, 10~^, 100) and obtain an output zq. Set xq = 

2) General Step {k = 1,2,...): 

a) Let D be the set of indices from 5'fc_i\Ji corresponding to the M components of Xfc_i with the 
smallest absolute value. Let E be the set of indices from S^_^ n J2 corresponding to the P components 
of V/(xfc_i) with the highest absolute value. 

b) Set 5 = Sk-i- For each i ^ D and for each j ^ E make a swap between the indices 

5 = (5fc_i\{i})U{i}. 

Invoke DGN with input (Aj, Cj, S*, 10^^, 100) and obtain an output z. Set x = U5Z. Advance T: 

r ^ r + 1. 

If /(x) < /(xfc-i), then set Sk = S,:Kk = x, advance k and goto 2. a. 

c) If none of the swaps resulted with a better objective function value, then STOP. The output is x = Xfe_i 
and T. 



random indices chosen to satisfy the support constraints Ji Q S Q J2- Then, at each iteration a swap between 
a support and an off-support index is performed such that the resulting solution via the DGN method improves 
the objective function. Since at each iteration only two elements are changed (one in the support and one in the 
off-support), this is a so-called "2-opt" method (see IITBI ). The swaps are always chosen to be between the M 
support indices corresponding to components in the current iterate ^k-i with the smallest absolute value and the 
off-support indices corresponding to the P components of V/(xfc_i) = 4 ^^(x^_j^AjXfc_i — Cj)AjXfc_i with 
the largest absolute value. This process continues as long as the objective function decreases and stops when no 
improvement can be made. 

A detailed description of the method is given in Algorithm [2] Note that the maximum number of swaps at each 
iteration of the method is MP. In our implementation the parameters were chosen to be M = 4, P = 8, which 
proved empirically to work well. We also note that the order of the elements in the indices set D and E was chosen 
to be random. 

C. The GESPAR Algorithm 

The 2-opt method can have the tendency to get stuck at local optima points. Therefore, our final algorithm, 
which we call GESPAR, is a restarted version of 2-opt. The 2-opt method is repeatedly invoked with different 
initial random support sets until the resulting objective function value is smaller than a certain threshold (success) 
or the number of maximum allowed total number of swaps was passed (failure). A detailed description of the 
method is given in Algorithm |3] One element of our specific implementation that is not described in Algorithm [3] 
is the incorporation of random weights added to the objective function, giving randomly different weights to the 
different measurements. Namely, the objective function used is actually chosen as /(x) = J2iLi ' (^^AjX — Cj)^ 
with Wj = 1 or 2 with equal probability. The random generation of weights is done each time the DGN procedure 
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Algorithm 3 GESPAR 
Input: (Ai,Ci,T,ITER). 

Aj G R^^^ , i = 1,2, . . . , N - symmetric matrices. 
CieR,i = l,2,...,N. 
T - threshold parameter. 

ITER - Maximum allowed total number of swaps. 
Output: X - an optimal (or suboptimal) solution of ([3]). 
Initialization. Set C = 0, /c = 0. 

• Repeat 

Invoke the 2-opt method with input (Aj, Cj, 4, 8) and obtain an output x and T. Set x^ = x, C = C + T and 

advance k: k -i^ k + 1. 

Until /(x) < T or C > ITER. 

• The output is x^ where I = argmin /(x^). 

m=0,l...fc-l 



is invoked. We observed that this modification reduced the probability of the 2-opt procedure to get stuck in 
non-optimal points. 

IV. Numerical Simulations 

In order to demonstrate the performance of GESRA^R, we conduct several numerical simulations. The algorithm is 
compared to other existing methods, and is evaluated in terms of signal-recovery accuracy, computational efficiency, 
and robustness to noise. 

A. Signal-recovery Accuracy 

In this subsection we examine the accuracy of GESR\R as a function of sparsity - i.e. the recovery success rate 
as a function of the number of non-zero elements in the signal. A runtime comparison of the compared methods 
is also performed. 

1 ) Simulation details: We choose x as a random vector of length n. The vector contains uniformly distributed 
values in the range [—4, —3] U [3, 4] in s randomly chosen elements. The point DFT of the signal is calculated, 
and its magnitude is taken as y, the vector of measurements. The 2n — 1 point correlation is also calculated. In 
order to recover the unknown vector x, the GESR\R algorithm is used with r = 10^^ and T = 20000, as well 
as two other algorithms for comparison purposes: An SDP based algorithm (Algorithm 2, |9|.), and an iterative 
Fienup algorithm with a sparsity constraint | |12)| . In our simulation n = 64 and N = 128. The Sparse-Fienup 
algorithm is run using 100 random initial points, out of which the chosen solution is the that best matches the 
measurements. Namely,x is selected as the s sparse output of the Sparse-Fienup algorithm with the minimal cost 
/(x) = E£i(|Fix|2 - y2)2 ojjt of thelOO runs. 

2) Simulation Results: Signal recovery results of the numerical simulation are shown in Fig. [T] where the 
probability for successful recovery is plotted for different sparsity levels. Successful recovery probability is defined 
as the ratio of correctly recovered signals x out of 100 signal-simulations. In each simulation both the support and 
the signal values are randomly selected. The three algorithms (GESPAR, SDP and Sparse-Fienup) are compared. 
The results clearly show that GESPAR outperforms the other algorithms in terms of probability of successful 
recovery - over 90% successful recovery up to s = 15, vs. s = 8 and s = 7 in the other two algorithms. 

Average runtime comparison of the three algorithms is shown in Table |l] The runtime is averaged over all 
successful recoveries. In all three algorithms the same vector sizes were used - i.e. n = 64, = 128. The 
computer used to solve all problems has an Intel 15 CPU and 4GB of RAM. As seen in the table, the SDP based 
algorithm is significantly slower than the other two, and the Fienup based algorithm, although fast, leads to a much 
lower success rate than GESPAR, which is fast and accurate. 
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TABLE I 
Runtime comparison 









SDP 


1 Sparse-Fienup 


GESPAR 






recovery 


runtime 


1 recovery | 


runtime 


recovery 


runtime 


s 


= 3 


0.93 


1.32 sec 


1 0.98 1 


0.09 sec 


1 


0.12 sec 


s 


= 5 


0.86 


1.78 sec 


1 0.97 1 


0.12 sec 


1 


0.12 sec 


s 


= 8 


0.9 


3.85 sec 


1 0.82 1 


0.50 sec 


1 


0.23 sec 



B. Sensitivity to exact sparsity knowledge 

Since the exact value of the signal's sparsity s is sometimes not known, the performance of GESPAR is examined 
when only an upper limit on s is known. 

1 ) Simulation Details: A GESPAR simulation is run twice: Once with s known exactly at each realization, and 
once with only an upper limit on s being known. The upper limit is taken as 32. The other simulation settings are 



the same as in Section |IV-AH this time with T = 5000. 

2 ) Simulation Results: Figure [2] shows the probability for successful recovery of the two simulations. The rather 
loose upper limit on s does not seem to affect the results significantly. 

C. Robustness to noise 

The performance of GESPAR as a function of SNR is examined, and compared with an existing method - sparse 
Fienup lfT2l . The SDP based method that is presented in ||9l is not designed to deal with noise as presented, so 
it was left out of the current comparison. The SDP method presented in flOl is designed to deal with random 
measurements and does not produce comparable results from direct Fourier measurements, and therefore is not 
compared with in this work. 



1) Simulation Details: As in Section IV- A 1 we choose x as a vector of length n, with s randomly chosen 
elements containing uniformly distributed values, and the N point Fourier magnitude, y, is calculated. Now, white- 
gaussian noise v is added to the measurements, at different SNR values, defined as: SNR = 201ogj|^. In order 
to recover the unknown vector x, the GESPAR algorithm is used with r = 10"'^ and T = 10000, as well as the 
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sparse -Fienup algorithm, for comparison purposes. In our simulation n = 64 and N = 128. The sparse-Fienup 
algorithm is run with a maximum of 1000 iterations, and with 100 random initial points. 

Note that even with little noise, the information on the support obtained by the zeros of the autocorrelation 
is no longer available. This is due to the fact that in the presence of noise, there will be no true zeros in the 
measured (or calculated) autocorrelation. In this case, one might try to threshold the autocorrelation values, rendering 
small autocorrelation values as zeros. However, this might result in zeroing of small (yet non-zero) values of 
the true autocorrelation function. Therefore, in the noisy case, we do not use support information obtained by 
the autocorrelation function in GESPAR, namely J2 = {1,2, .. . ,n}. The noiseless case was simulated with the 
autocorrelation support information. 

2) Simulation Results: Figure [3] shows the normalized mean squared reconstruction error (NMSE), defined as 
NMSE = as a function of sparsity, for different SNR values. Each point represents an average over 300 

different random realizations. The performance under different SNR values is plotted for GESPAR (full lines), 
and for sparse-Fienup (dashed-lines). The performance of GESPAR naturally improves as SNR increases, and it 
outperforms the sparse-Fienup in terms of noise-robustness. 



D. Scalability 

As the main advantage of GESPAR over SDP based methods is its ability to solve large problems efficiently, its 
performance for different vector sizes is examined. 

1) Simulation Details: GESPAR is simulated for various values of n (64,128,256,512). In all simulated cases 



N = 2n. The other simulation parameters are as in Section IV-Bl 



2 ) Simulation Results: The recovery probability vs. sparsity s for different vector lengths is shown in Fig. |4] 
The maximal sparsity allowing succesful recovery is shown to increase with vector length n. The maximal value 
of s that is recoverable with high probability seems to scale like n^/^, which is consistent with the same scaling 
observation presented in ||9l. The mean reconstruction time for a signal with n = 512, s = 35 from N = 1024 
measurements, allowing T = 5000 replacements, is 25.92 seconds. 
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Fig. 3. Normalized MSB vs. sparsity level. The performance is plotted for several SNR values for GESPAR (full lines) and Sparse-Fienup 
(dashed lines). 
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Fig. 4. Scalability - recovery probability as a function of signal sparsity, for various vector lengths (64,128,256,512) 

V. Conclusion 

We proposed and demonstrated GESPAR - a fast algorithm to recover a sparse vector from its Fourier magnitude. 
We showed via simulations that GESPAR outperforms alternative approaches suggested for this problem. The 
algorithm does not require matrix-lifting, and therefore is potentially suitable for large scale problems such as 2D 
images. 
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